Prévision de la performance énergétique des bâtiments¶
La performance énergétique est un enjeu majeur de la réduction de la consommation d’énergie et des émissions de gaz à effet de serre. La menace du changement climatique pèse sur l’humanité et la transition énergétique vise notamment à améliorer cette performance énergétique des bâtiments. Depuis plus d’une décennie, le diagnostic de performance énergétique vient faciliter ce travail en apportant des données clés.
L'étude menée se distingue en 2 phases :
- Exploration statistique pour comprendre les différentes variables explicatives du jeu de données et observer des éventuelles corrélations.
Modélisation des performances énergétiques de nouveaux batiments à l'aide d'algorithmes prédictifs (Random Forest, SVM, etc). Il est important de noter que nous allons traiter le problème selon 2 aspects : Classification et Régression.
Importation des packages¶
library(corrplot)
library(ggplot2)
library(lattice)
library(rpart)
library(partykit)
library(randomForest)
library(ROCR)
library(missForest)
library(reshape2)
library(FactoMineR)
library(factoextra)
#library(rgl)
library(plotly)
library(dplyr)
library(gridExtra)
library(bestglm)
library(MASS)
library(VGAM)
library(e1071)
library(glmnet)
library(broom)
library(gbm)
library(bestglm)
library(adabag)
library(caret)
Importation des données¶
data.complet <- read.csv('DataEnergy.csv')
x <- data.complet
head(x)
Descriptif des variables¶
- Relative compactness : compacité relative (propriété liée à la forme du batiment et sa surface).
- Surface area : surface totale du batiment.
- Wall area : surface des murs.
- Roof area : surface du toit.
- Overal height : hauteur du batiment (variable qualitative ordonnée: 3,5m et 7m).
- Orientation : orientation du batiment (North, East, South, West).
- Glazing area : surface totale des vitrages du batiment.
- Glazing area distr : orientation des vitrages du batiment (55% East, 55% South, etc).
- Energy : quantifie les performances énergétique du batiment.
- Energy efficiency : lettre de classification des performances énergétique du batiment.
Préparation des variables¶
str(x)
# Transformation des variables qualitatives en variables catégorielles
x$Glazing.area.distr <- as.factor(x$Glazing.area.distr)
levels(x$Glazing.area.distr) <- c("No Glazing", "55% North", "55% East",
"55% South", "55% West", "Uniform")
x$Energy.efficiency <- factor(x$Energy.efficiency, ordered=TRUE)
summary(x)
On remarque des valeurs négatives pour la variable Glazing Area, qui représente pourtant la surface totale des vitrages. En fait, ces valeurs négatives sont dues à un bruit ajouté. Regardons comment sont réparties ces valeurs négatives afin de leur appliquer une correction.
summary(x[which(x$Glazing.area<0),]$Glazing.area.distr)
D'après le resultat ci-dessus, toutes les valeurs de Glazing.area inférieures à 0 appartiennent à la catégorie No glazing. On va donc forcer leur valeur à 0.
x$Glazing.area[x$Glazing.area < 0] <- 0
On remarque aussi que la variable Overall.height prend seulement deux niveaux ordonnés : 3,5m et 7m. On va donc considérer cette variable comme catégorielle ordonnée.
x$Overall.height <- as.ordered(x$Overall.height)
A. Analyse des données¶
1. Statistiques descriptives univariées¶
options(repr.plot.width = 4, repr.plot.height = 4)
histogram(x$Energy.efficiency, freq=FALSE, xlab="Energy Efficiency",
main="Proportion du nombre de batiments \npar classe énergétique")
La répartition des données par classe énergie est plutot uniforme. On observe cependant que la classe A est majoritaire dans ce jeu de données.
options(repr.plot.width = 6, repr.plot.height = 6)
boxplot(x[,-c(5,6,8,10,11)], col='lightblue', xaxt = "n",
main="Box plot des variables quantitatives")
grid()
text(seq_along(x[,-c(5,6,8,10,11)]), par("usr")[3] - 0.5, labels = names(x[,-c(5,6,8,10,11)]), srt = 45, adj = 1, xpd = TRUE);
Les variables ne sont pas toutes de même ordre de grandeur mais les distributions semblent homogènes. Etant donné les ordres de grandeurs bien différents, une attention particulière sera accordée lors de la réalisation de l'ACP.
Il faut maintenant étudier plus précisemment la distribution de ces variables.
par(mfrow=c(2,3))
options(repr.plot.width = 9, repr.plot.height = 6)
hist(x$Glazing.area, xlab='Glazing Area', main='Glazing Area Distribution',col='lightblue')
hist(x$Relative.compactness,xlab= 'Relative Compactness' ,main = 'Relative Compactness Distribution',col='lightblue')
hist(x$Roof.area,xlab= 'Roof Area' ,main = 'Roof Area Distribution',col='lightblue')
hist(x$Surface.area,xlab= 'Surface Area' ,main = 'Surface Distribution',col='lightblue')
hist(x$Wall.area,xlab= 'Wall Area' ,main = 'Wall Area Distribution',col='lightblue')
hist(x$Energy,xlab= 'Energy' ,main = 'Energy Distribution',col='lightblue')
Il n'y a pas de transformation de variables nécessaires. On remarque tout de même que la variable Roof Area est inéquitablement répartie.
options(repr.plot.width = 4, repr.plot.height = 4)
pie(table(x$orientation), main = "Proportion du nombre de batiment \n selon l'orientation")
La répartition des orientations des bâtiments est uniforme pour le jeu de données.
options(repr.plot.width = 6, repr.plot.height = 3)
histogram(x$Glazing.area.distr,xlab= 'Glazing area distr' ,
main = 'Proportion du nombre de batiments \n par orientation des vitrages',col='gold')
La répartition des orientations des vitrages est équitable entre les différentes classes (55% Nord, 55% Est, 55% Sud, 55% Ouest et sans vitrage). Cependant, les batiments sans vitres sont minoritaires.
2. Statistiques descriptives multivariées¶
options(repr.plot.width = 6, repr.plot.height = 5)
M <- cor(x[,-c(5,6,8,10,11)])
corrplot(M)
Relative compactness est très fortement correlé négativement avec Surface area (-0.98) et Roof Area (-0.87). Il doit exister un lien numérique entre ces variables que nous étudierons par la suite.
Roof Area est aussi fortement corrélé avec Energy.
Il est difficile de tirer des conclusions avec le graphique des corrélations. Regardons la matrice de scatterplots.
options(repr.plot.width = 6, repr.plot.height = 6)
pairs(x[,-c(5,6,8,10,11)])
On remarque une relation linéaire entre Relative.compactness et Surface.area.
On voit que Glazing.area, qui est une variable quantitative, est finalement répartie selon 4 intervalles.
defaultW <- getOption("warn")
options(warn = -1)
fig <- plot_ly(x, x = ~Roof.area, y = ~Wall.area, z = ~Surface.area, color = ~Energy.efficiency)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Roof Area'),
yaxis = list(title = 'Wall Area'),
zaxis = list(title = 'Surface Area')))
fig
options(warn = defaultW)
## On vérifie la relation Surface.area = Wall.area + 2*Roof.area
a=x$Surface.area-(x$Wall.area+2*x$Roof.area)
cat("Nombre de valeurs pour lesquelles la relation 'Surface.area = Wall.area + 2*Roof.area' n'est pas vérifié : ", length(which(round(a,3)!=0)))
La surface est calculée en sommant la surface des murs et en ajoutant deux fois la surface du toit.
La relation étant vérifiée pour tous les points de notre jeu de donnée, nous allons pouvoir supprimer une ou plusieurs variables étant donné que les informations sont redondantes.
options(repr.plot.width = 18, repr.plot.height = 7)
par(mfrow=c(1,3))
boxplot(Wall.area ~Energy.efficiency, data=x, col='lightblue', main="Wall Area vs Energy Efficiency")
boxplot(Roof.area ~Energy.efficiency, data=x, col='lightgreen', main="Roof Area vs Energy Efficiency")
boxplot(Surface.area ~Energy.efficiency, data=x, col='lightyellow', main="Surface Area vs Energy Efficiency")
On remarque sur ces graphiques que les classes A, B et C ont une surface de toît élevée. De la même façon, ces classes ont une surface totale élevée.
La tendance est inversée pour les classes D, E, F et G.
Pour continuer l'analyse, on va donc supprimer la variable Roof Area, d'une part car son information est contenue dans la variable Surface Area et d'autre part car sa distribution est inéquitablement répartie.
On retire donc la variable Roof Area.
x <- dplyr::select(x,-Roof.area)
options(repr.plot.width = 10, repr.plot.height = 5)
par(mfrow=c(1,2))
boxplot(x$Overall.height~ x$Energy.efficiency, main ="Hauteur \n selon la classe Energetique", xlab="Classe Energetique", ylab="Hauteur",col='red')
boxplot(x$Relative.compactness~ x$Energy.efficiency, main ="Relative compactness \n selon la classe energetique", xlab="Classe Energetique", ylab="Hauteur",col='darkred')
On remarque sur le graphique que les classes A, B et C ont majoritairement une hauteur de 3m50 tandis que les classes D, E, F et G ont principalement une hauteur plus élevée (7m).
De la même façon, on remarque que la valeur de Relative Compactness est plus faible pour les bonnes classes énergetiques par rapport aux moins bonnes.
options(repr.plot.width = 12, repr.plot.height = 5)
h1 = histogram(~Energy.efficiency | orientation , data=x, xlab= "Energy Efficiency",main="Orientation vs Class Energy",)
h2 = histogram(~Energy.efficiency | Glazing.area.distr, data=x,main="Glazing area distr vs Class Energy", xlab= "Energy Efficiency")
grid.arrange(h1, h2, nrow=1, ncol=2)
#histogram(~orientation | Energy.efficiency , data=x, xlab= "Energy Efficiency",main="Orientation vs Class Energy",)
#mosaicplot(table2, main="Glazing area distr vs Class Energy")
Les classes d'énergie sont uniformement réparties selon l'orientation du bâtiment et selon l'orientation des vitrages. L'orientation des bâtiments n'influe donc probablement pas l'appartenance à une certaine classe d'énergie.
En revanche, lorsqu'il y a pas de vitrage, la classe A est clairement majoritaire. On peut faire le lien avec la réalité où on suppose qu'un batiment sans vitrage est mieux isolé.
3. Analyse en Composantes Principales¶
On enlève la variable Energy car c'est celle que l'on veut prédire. On garde seulement les variables quantitatives.
Remarque : Les variables n'étant pas dans les mêmes unités, nous faisons une ACP centrée réduite.
X_acp <- dplyr::select(x,c(Relative.compactness, Surface.area, Wall.area, Glazing.area,Energy.efficiency))
res.pca <- PCA(X_acp, scale.unit=TRUE, quali.sup=5, graph = FALSE)
#fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50), new.plot=FALSE, graph.type="classic")
lbls = paste("Dimension",1:4, ":", round(res.pca$eig[,2],2), "%" )
coul = c('palegreen','salmon','lightyellow','lightblue','black')
options(repr.plot.width = 8, repr.plot.height = 9)
par(mfrow=c(2,1))
barplot(res.pca$eig[,2],names.arg=paste("Dim.",1:nrow(res.pca$eig)), xlab="Dimensions", ylab = "% of Variance", col = coul,main="Pourcentage de variance expliquée par dimension")
legend(x='topright',legend=lbls, fill= coul, bty='n', cex=0.8, text.font=0.5)
boxplot(res.pca$ind$coord, outlier=TRUE, col=coul, main = "Distribution des premières composantes")
Au vu du pourcentage de variance expliqué et des ditributions des trois premières dimensions, nous sommes tentés de les garder.
options(repr.plot.width = 15, repr.plot.height = 6)
plot1 <- plot(res.pca,choix="ind", label="none", new.plot=FALSE, graph.type="ggplot",habillage=5)
plot2 <- plot(res.pca,choix="var", new.plot=FALSE, graph.type="ggplot")
grid.arrange(plot1, plot2, nrow=1, ncol=2)
On observe des clusters dans le graphe des individus qui ne reflètent pas exactement les classes. Cependant on peut séparer les "bonnes classes" (A, B et C) des moins bonnes classes (D,E,F et G) par la première composante principale. En effet, les "bonnes classes" semblent se distinguer des "mauvaises classes", par une valeur positive sur l'axe 1. Ce dernier étant porté par Surface Area positivement et par Relative Compactness négativement.
La dimension 2 est principalement portée par Glazing Area. On ne distingue pas les différentes classes grâce à cette composante (et donc à cette variable).
options(repr.plot.width = 15, repr.plot.height = 6)
plot3 <- plot(res.pca,axes=c(2,3),choix="ind", label="none", new.plot=FALSE, graph.type="ggplot",habillage=5)
plot4 <- plot(res.pca,axes=c(2,3),choix="var", new.plot=FALSE, graph.type="ggplot")
grid.arrange(plot3, plot4, nrow=1, ncol=2)
On remarque encore ici des clusters selon l'axe 2; ces derniers sont sans doute dus à la variable Glazing Area qui est disparate comme nous l'avons souligné précedemment.
On distingue tout de même que les bonnes classes ont tendance à avoir une valeur négative sur l'axe 3, porté par Wall Area. La variable Wall Area permet sans doute elle aussi de distinguer certaines classes, plus Wall Area est faible, plus la classe Energétique est bonne.
Les conclusions de l'ACP correspondent aux caractéristiques montrées précedemment sur l'impact de certaines variables sur la consomation Energetique.
4. Clustering¶
x_num=dplyr::select(x, Relative.compactness, Wall.area,Surface.area, Glazing.area)
x_num_cr=scale(x_num)
stacked_cluster = hclust(dist(x_num_cr), method = "ward.D")
plot(stacked_cluster)
# Elbow method
fviz_nbclust(x_num_cr, kmeans, method = "wss") +
geom_vline(xintercept = 2, linetype = 2)+
labs(subtitle = "Elbow method")
Ces deux méthodes nous suggèrent 2 clusters, or on aimerait qu'il y en ai 7 (correspondant aux 7 classes énergétiques)
set.seed(10)
km.out2 = kmeans(x_num_cr,centers=2)
km.out7 = kmeans(x_num_cr,centers=7)
plot.kmean2<-fviz_cluster(km.out2, x_num_cr, ellipse.type = "norm", main='Kmeans with 2 clusters')
plot.kmean7<-fviz_cluster(km.out7, x_num_cr, ellipse.type = "norm", main='Kmeans with 7 clusters')
grid.arrange(plot.kmean2, plot.kmean7, nrow=1, ncol=2)
cat('Table de contingence avec 2 clusters :')
print(table(km.out2$cluster,x$Energy.efficiency))
cat(' \n\n ')
cat('Table de contingence avec 7 clusters :')
print(table(km.out7$cluster,x$Energy.efficiency))
Kmeans ne semble pas adapté à notre problème. Cependant, cet algorithme suggère de regrouper les classes A,B et C ensemble ainsi que D,E,F,G (bonne consommation et mauvaise consommation). Nous nous intéresserons sur cela lors d'une partie Bonus.
B. Séparation des données :¶
On sépare les données en un ensemble d'apprentissage (train) et un ensemble test. Le modèle est entrainé sur l'échantillon train et il est évalué sur l'échantillon test.
Cette étape est nécessaire pour permettre d'évaluer le modèle déterminé avec les échantillons d'entrainement. En effet, on peut par exemple détecter un overfitting grâce à l'échantillon test (un très bon score avec l'échantillon d'entrainement ne signifie pas toujours que c'est un très bon modèle).
## 75% of the sample size
x_fin<-x
train_ratio=0.7
smp_size <- floor(train_ratio * nrow(x))
## set the seed to make your partition reproducible
set.seed(2512)
train_ind <- sample(seq_len(nrow(x)), size = smp_size)
train <- x[train_ind, ]
test <- x[-train_ind, ]
cat("Nombre d'individus dans l'échantillon train : ", nrow(train), '\n')
cat("Nombre d'individus dans l'échantillon test : ", nrow(test))
options(repr.plot.width = 5, repr.plot.height = 4)
histogram(train$Energy.efficiency, freq=FALSE, xlab="Energy Efficiency",
main="Proportion du nombre de batiments \npar classe énergétique dans le train set")
On remarque une plus forte proportion d'individus de la classe A. Nos modèles auront donc tendance à être meilleurs pour la classe A.
Nous avons choisi des métriques de performance qui tiennent compte de ce fait.
C. Modélisation¶
Définition des métriques¶
Pour évaluer les performances d'un modèle de classification, il faut définir des mesures traduisant son efficacité:
Quatres scores ont été utilisées pour comparer les différents modèles :
Score accuracy : Quotient entre les bien prédits et le nombre total de prédictions (fonction perf_model). Le défaut de cette première métrique est qu'elle ne prend pas en compte l'écart de classe quand on se trompe. En effet, imaginons que la vraie valeur soit la classe A ; l'erreur est bien plus grave si on prédit une consommation de type G que si l'on prédit une consommation de type B.
Root Mean Squared Error (RMSE) : Uniquement utilisé dans le cadre de la régression
Mean Absolute Error (MAE) : Cette métrique est obtenue en utilisant la formule trouvée dans le papier Measuring the performance of ordinal classification (Aime S. Cardoso & Ricardo Sousa). Bien qu'associé en temps normal à des regréssions, cette mesure est ici adaptée pour prendre en compte les écart de classe lors d'une classification. (fonction MAE()). Plus le score est faible, meilleur est la classification.
Macro F1 calcule la métrique F1 pour chaque classe individuellement (comme pour une classification binaire où la valeur TRUE correspondrait à la classe que l'on souhaite prédire la valeur FALSE regrouperait toutes les autres classes.). On fait ensuite la moyenne de toutes ces mesures F1. Il est important de noter que cette métrique ne tient pas en compte des classes non uniformément réparties.
NB : Nous gardons aussi toutes les valeurs individuelles de mesures F1 et nous les afficherons à la fin pour comparer nos modèles. Certaines méthodes (linéaire ou non-linéaires) seront en effet sans doutes meilleures pour prédire certaines classes que d'autres.
Toutes les métriques pour la multiclassification peuvent être trouvées ici :
http://www.inescporto.pt/~jsc/publications/journals/2011JaimeIJPRAI.pdf
## Définition des scores
perf_model <- function(table_contingence){
return( sum(diag(table_contingence))/sum(table_contingence))
}
MAE <- function(table){
lignes = nrow(table)
cols = ncol(table)
N = sum(table)
MSE=0
for (i in seq(1,lignes)){
for (j in seq(1,cols)){
MSE = MSE + (table[i,j]*abs(i-j))
}
}
return(MSE/N)
}
macro_F1 <- function(table){
lignes = nrow(table)
cols = ncol(table)
diag = diag(table)
rowsums = apply(table, 1, sum)
colsums = apply(table, 2, sum)
precision = diag / colsums
recall = diag / rowsums
f1 = 2 * precision * recall / (precision + recall)
return(f1)
}
#Definition des tableaux de score :
#Ces fonctions permettent de regrouper dans un tableau les différents scores pour chaque modèle
Tab_score_class = as.data.frame(setNames(replicate(4,numeric(0), simplify = F),
c("Modele","Accuracy","MAE", "Macro F1") ))
Tab_score_reg = as.data.frame(setNames(replicate(5,numeric(0), simplify = F),
c("Modele","Accuracy","MAE", "Macro F1","RMSE") ))
Tab_F1_class = as.data.frame(setNames(replicate(8,numeric(0), simplify = F),
c("Modele","A","B","C","D","E","F","G") ))
Tab_F1_reg = as.data.frame(setNames(replicate(8,numeric(0), simplify = F),
c("Modele","A","B","C","D","E","F","G") ))
#pred est la liste des prédictions sur le jeu de test
#quali est le type de la prediction
Compute_Error <- function(pred,quali=TRUE, name_model=""){
if (quali){
table_result = table(pred.reg = pred , observations = test$Energy.efficiency)
cat("Table de contingence de", name_model, ": \n")
print(table_result)
score1 = perf_model(table_result)
score2 = MAE(table_result)
score3 = mean( macro_F1(table_result) )
#on cree une df avec les mêmes noms de colonnes
to_add <- data.frame(Modele=name_model, Accuracy=score1,MAE=score2, MacroF1=score3)
Tab_score_class = rbind(Tab_score_class,to_add) # on concatène ensuite
return(Tab_score_class)
}
else{
pred_class <- classify(pred)
table_result = table(pred.reg = pred_class , observations = test$Energy.efficiency)
cat("Table de contingence de", name_model, ": \n")
print(table_result)
score1 = perf_model(table_result)
score2 = MAE(table_result)
score3 = mean( macro_F1(table_result) )
rmse = sqrt( sum( (pred - test$Energy)^2 ) /nrow(test) )
#on cree une df avec les mêmes noms de colonnes
to_add <- data.frame(Modele=name_model, Accuracy=score1,MAE=score2, MacroF1=score3, RMSE=rmse)
Tab_score_reg = rbind(Tab_score_reg,to_add) # on concatène ensuite
return(Tab_score_reg)
}
}
Compute_F1 <- function(predict,quali=TRUE, name_model=""){
if (quali){
table_result = table(pred.reg = predict , observations = test$Energy.efficiency)
temp = macro_F1(table_result)
add_F1= data.frame(Modele = name_model,"A"= temp[[1]],"B"= temp[[2]], "C"= temp[[3]], "D"=temp[[4]],
"E"=temp[[5]],"F"=temp[[6]],"G"=temp[[7]])
Tab_F1_class = rbind(Tab_F1_class,add_F1)
return(Tab_F1_class)
}
else{
pred_class <- classify(predict)
table_result = table(pred.reg = pred_class , observations = test$Energy.efficiency)
temp = macro_F1(table_result)
add_F1= data.frame(Modele=name_model,"A"=temp[1],"B"=temp[2], "C"=temp[3], "D"=temp[4],
"E"=temp[5],"F"=temp[6],"G"=temp[7])
Tab_F1_reg = rbind(Tab_F1_reg,add_F1)
return(Tab_F1_reg)
}
}
1. Classification¶
Dans cette partie, nous nous concentrons uniquement sur le problème de classification, qui consiste à prédire la classe énergétique d'un batiment. La variable quantitative Energy est donc supprimée de cette partie.
train_class <- dplyr::select(train,-Energy)
test_class <- dplyr::select(test, -Energy)
x_train_class<-dplyr::select(train_class, -Energy.efficiency)
y_train_class<- train_class$Energy.efficiency
x_test_class<-dplyr::select(test_class,-Energy.efficiency)
y_test_class<-test_class$Energy.efficiency
a. Régression logistique polythomiale :¶
Nous allons réaliser une régression logistique polynomial ordonnée aussi appelée "polythomial regression".
Il y a deux moyens de le réaliser : additive logits and adjacents logits.
Mise en place du modèle:
#Additive Simple:
vglm.c <- vglm(Energy.efficiency ~., data=train_class, family=cumulative(parallel=T, reverse=F))
Résultats :
p <- predict(vglm.c, newdata = test_class, type = "response") #matrice de proba d'appartenance à chaque classe
vglm.c.pred <- apply(p,1, which.max)
Tab_score_class <-Compute_Error(vglm.c.pred,quali=TRUE,name_model = "Naive Logit");Tab_score_class
Tab_F1_class <- Compute_F1(vglm.c.pred,quali=TRUE,name_model = "Naive Logit")
Commentaire :
Le score de classification est relativement faible, mais il n'est pas mauvais. On remarque par ailleurs que les erreurs de classification sont majoritairement dans les classes voisines. Ce résultat n'est pas choquant et c'est un bon début. Dans la partie suivante, nous allons essayer de pénaliser cette régression.
b. Regression logistique pénalisée :¶
Mise en place du modèle :
x.mat.train <- model.matrix(Energy.efficiency ~ . -1 , data = train_class) #permet de gérer les variables catégorielles
lasso.c.cv <- cv.glmnet(y = y_train_class, x = x.mat.train,family="multinomial")
plot(lasso.c.cv)
paste("Estimation CV de lambda : ", round(lasso.c.cv$lambda.min,3))
Résultats :
x.mat.test <- model.matrix(Energy.efficiency ~ . -1, data = test_class)
lasso.c.cv.pred <- predict(lasso.c.cv,s="lambda.min", newx = x.mat.test,type="class") #valeur predite sur le jeu de train
Tab_score_class <-Compute_Error(lasso.c.cv.pred,quali=TRUE,name_model = "Penalised Logit");Tab_score_class
Tab_F1_class <- Compute_F1(lasso.c.cv.pred,quali=TRUE,name_model = "Penalised Logit")
Commentaires :
La pénalisation LASSO permet de diminuer le nombre de variables utilisées tout en dégradant très peu l'information. Dans ce modèle une contrainte est ajoutée sur les paramètres. La pénalisation LASSO a donc une moins bonne accuracy que le modèle linéaire car le modèle linéaire minimise uniquement les moindres carrés (et ne pénalise pas les variables).
Nous avons donc un moins bon modèle mais ce modèle est plus parcimonieux : il contient moins de variables.
Essayons maintenant des méthodes non-linéaires.
c. Arbre de décision¶
Mise en place du modèle optimisé
options(repr.plot.width = 15, repr.plot.height = 10)
tree.c=rpart(Energy.efficiency~.,data=train_class,control=rpart.control(cp=0.001))
xmat=xpred.rpart(tree.c, xval = 10)
#levels(datappr[,"Survived"]) <- c(1,2)
xerr <- (xmat-as.numeric(y_train_class))^2
CVerr=apply(xerr,2,sum)
cp_opti = as.numeric(attributes(which.min(CVerr))$names);
cat('cp opti :', cp_opti)
tree.c.opti=rpart(Energy.efficiency~.,data=train_class,control=rpart.control(cp=cp_opti))
options(repr.plot.width = 15, repr.plot.height = 10)
plot(tree.c.opti)
text(tree.c.opti)
Commentaire : On remarque que la variable Overall.height permet de classifier les extrêmes (classe A,B vs D,E,F,G), la délimitation dela classe C n'étant pas aussi nette.
Résultats :
tree.c.opti.reg <- predict(tree.c.opti, newdata = test_class, type = "class")
Tab_score_class <-Compute_Error(tree.c.opti.reg,quali=TRUE,name_model = "Tree opti");Tab_score_class
Tab_F1_class <- Compute_F1(tree.c.opti.reg,quali=TRUE,name_model = "Tree Opti")
Commentaire :
L'apport de cette méthode est significative car notre prédiction s'améliore pour chacune des métriques étudiées. L'utilisation de méthodes non-linéaires semble pour l'instant être une bonne piste. Nous allons donc par la suite en étudier d'autres. On se dirige donc vers un RandomForest qui devrait prodiguer d'encore meilleurs résultats.
d. Random Forest¶
Mise en place du modèle :
rf.c<-randomForest::randomForest(x = x_train_class, y = y_train_class,
xtest = x_test_class, ytest = y_test_class,
importance = TRUE)
Résultats :
rf.c.pred <- rf.c$test$predicted
Tab_score_class <-Compute_Error(rf.c.pred,quali=TRUE,name_model = "Random Forest");Tab_score_class
Tab_F1_class <- Compute_F1(rf.c.pred,quali=TRUE,name_model = "RF");
Commentaire : Nous remarquons dans un premier temps que le taux de classification est meilleur que celui obtenu avec l'arbre de décision.
Pour améliorer le taux de classification de Random Forest, nous pouvons essayer de modifier le paramètre mtry correspondant au nombre de variables testées à chaque division.
Utilisons la librairie caret.
Mise en place du modèle optimisé :
cvControl <- trainControl(method = "cv", number = 10)
mtryTrials <- train(x_train_class, y_train_class, method = "rf", tuneLength = 6,
trControl = cvControl, trace = FALSE)
mtryTrials
rf.c.opti <- randomForest::randomForest(x = x_train_class, y = y_train_class,
xtest = x_test_class, ytest = y_test_class,
mtry=mtryTrials$bestTune$mtry,importance = TRUE)
Résultats du modèle optimisé :
rf.c.opti.pred <- rf.c.opti$test$predicted
Tab_score_class <-Compute_Error(rf.c.opti.pred,quali=TRUE,name_model = "Random Forest Opti");Tab_score_class
Tab_F1_class <- Compute_F1(rf.c.opti.pred,quali=TRUE,name_model = "RF_Opti");
Commentaire :
Le modèle Random Forest retourne jusqu'ici les meilleurs résultats. La MAE est faible et les métriques F1 semblent être élevées.
On voit donc que l'apport du modèle optimisé pour le nombre de variables retenues aléatoirement à chaque noeud (mtry) n'est pas forcément significatif et ajoute un cout numérique supplémentaire.
e. Boosting¶
Mise en place du modèle :
defaultW <- getOption("warn")
options(warn = -1)
gbm.c <- gbm(Energy.efficiency ~ ., data=train_class,
distribution = "multinomial",cv.folds = 10,shrinkage = .01,
n.minobsinnode = 10, n.trees = 200, verbose=FALSE)
options(warn = defaultW)
Résultats :
p = predict.gbm(object = gbm.c,
newdata = test_class,
n.trees = 200,
type = "response")
gbm.c.pred = colnames(p)[apply(p, 1, which.max)]
Tab_score_class <- Compute_Error(gbm.c.pred,quali = TRUE, name_model = "GBM"); Tab_score_class
Tab_F1_class <- Compute_F1(gbm.c.pred,quali=TRUE,name_model = "GBM");
Mise en place du modèle optimisé :
set.seed(2)
defaultW <- getOption("warn")
options(warn = -1)
cvControlRandom = trainControl(method = "cv", number = 10, search='random')
gbm.opti.c = train(x_train_class, y_train_class,method = "gbm", tuneLength = 8,
trControl = cvControl, verbose=FALSE)
options(warn = defaultW)
Résultats :
gbm.opti.c.pred = predict(gbm.opti.c, newdata = x_test_class)
Tab_score_class <- Compute_Error(gbm.opti.c.pred,quali = TRUE, name_model = "GBM opti"); Tab_score_class
Tab_F1_class <- Compute_F1(gbm.opti.c.pred,quali=TRUE,name_model = "GBM opti")
Commentaire :
Que ce soit le modèle optimisé ou non optimisé, utiliser la méthode de Gradient Boosting ne semble pas apporter de réelles améliorations par rapport à Random Forest, même avec des paramètres optimisés.
f. Support Vector Machine (SVM)¶
Mise en place du modèle :
svm.c=svm(Energy.efficiency~., data=train_class, kernel="radial")
Résultats :
svm.c.pred=predict(svm.c, newdata=test_class, type = "class")
Tab_score_class <- Compute_Error(svm.c.pred,quali = TRUE, name_model = "SVM"); Tab_score_class
Tab_F1_class <- Compute_F1(svm.c.pred,quali=TRUE,name_model = "SVM")
Mise en place du modèle optimisé :
defaultW <- getOption("warn")
options(warn = -1)
svm.c.tune <- tune.svm(Energy.efficiency~., data=train_class, type = "C",
kernel=c('polynomial','radial','linear'), degree=c(2,3,4),
cost=c(1,1.25,1.5,1.75,2),gamma = seq(0.1, 2, by = 0.2))
options(warn = defaultW)
(svm.c.tune$best.parameters)
Résultats du modèle optimisé :
svm.c.opti=predict(svm.c.tune$best.model, newdata=test_class, type = "class")
Tab_score_class <- Compute_Error(svm.c.opti, quali = TRUE, name_model = "SVM Opti");
Tab_score_class
Tab_F1_class <- Compute_F1(svm.c.opti,quali=TRUE,name_model = "SVM Opti")
Commentaire :
Tout comme GBM, la SVM n'apporte pas de meilleurs résultats de prédiction, la coût numérique de ces méthodes et de leur optimisation est d'ailleurs très élevé.
Bilan - Classification¶
cat("Bilan des scores de classification : \n " )
Tab_score_class
Pour l'interprétation, se référer à la partie Conclusion.
2. Régression¶
Dans cette partie, nous allons traîter le problème comme un problème de regression. En effet, la variable Energy est une variable quantitative continue et c'est à partir de cette variable que les classes sont formées :
- Les bâtiments ayant une consommation inférieure à 30 correspondent à une classe A
- Les bâtiments ayant une consommation comprise entre 30 et 35 correspondent à une classe B
- Les bâtiments ayant une consommation comprise entre 35 et 45 correspondent à une classe C
- Les bâtiments ayant une consommation comprise entre 45 et 55 correspondent à une classe D
- Les bâtiments ayant une consommation comprise entre 55 et 65 correspondent à une classe E
- Les bâtiments ayant une consommation comprise entre 65 et 75 correspondent à une classe F
- Les bâtiments ayant une consommation supérieure à 75 correspondent à une classe F
Ainsi nous allons prédire la variable continue Energy, puis attribuer une classe à la prédiction obtenue (selon la règle ci-dessus).
#New train and test sets : we remove Energy.efficiency
train_reg <- dplyr::select(train, -Energy.efficiency)
test_reg <- dplyr::select(test, -Energy.efficiency)
x_train_reg<-dplyr::select(train_reg, -Energy)
y_train_reg<-train_reg$Energy
x_test_reg<-dplyr::select(test_reg, -Energy)
y_test_reg<-test_reg$Energy
# permet de trouver la classe energetique à partir de la valeur numérique prédite
classify <- function(quantitative.predict){
predict.class = (quantitative.predict)
predict.class[quantitative.predict <=30] = 'A'
predict.class[quantitative.predict >30 & quantitative.predict <= 35] = 'B'
predict.class[quantitative.predict >35 & quantitative.predict <= 45] = 'C'
predict.class[quantitative.predict >45 & quantitative.predict <= 55] = 'D'
predict.class[quantitative.predict >55 & quantitative.predict <= 65] = 'E'
predict.class[quantitative.predict >65 & quantitative.predict <= 75] = 'F'
predict.class[quantitative.predict > 75] = 'G'
return(predict.class)
}
a. Regression linéaire¶
Mise en place du modèle :
glm.r <- glm(Energy ~., data=train_reg)
options(repr.plot.width = 8, repr.plot.height = 4)
par(mfrow=c(1,2))
plot(glm.r, which=1:2)
Commentaires : Le graphe précédent nous permet de visualiser si les hypothèses de linéarité et d'homoscédasticité du modèle.
Les résidus sont centrés mais il semble y avoir une dépendance entre la valeur prédite et l'erreur (quand valeur prédite élevée : erreur forte); graphe en forme de banane.
summary(glm.r)
Résultats :
glm.r.pred <- predict(glm.r, newdata=test_reg)
Tab_score_reg <- Compute_Error(glm.r.pred,quali = FALSE, name_model = "Regression Linéaire Naive");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.r.pred,quali=FALSE,name_model = "Reg Naive")
Commentaire :
Une régression linéaire naive a été réalisée : tous les paramètres ont été ajoutés au modèle.
La RMSE reste élevé pour distinguer toutes les classes. En effet, l'étendu des différentes classes est souvent égale à 10 (même parfois à 5 pour la classe B). La RMSE reste donc trop élevé lors de la séparation entre les différentes classes.
Les résultats obtenues avec la regression ne semble pas apporter une amélioration au modèle.
b. Régression linéaire pénalisée¶
Mise en place du modèle :
x.mat.train <- model.matrix(Energy ~ . - 1, data = train_reg) #permet de gérer les variables catégorielles
lasso.r <- glmnet(y = y_train_reg, x = x.mat.train)
options(repr.plot.width = 10, repr.plot.height = 7)
plot(lasso.r, xvar = "lambda", label = TRUE)
grid()
legend("bottomright",
legend = paste(1:ncol(x.mat.train), " - ", colnames(x.mat.train)), ncol=2, cex=0.75)
# choix du paramètre de régularisation par validation croisée
lasso.r.cv <- cv.glmnet(y = y_train_reg, x = x.mat.train)
plot(lasso.r.cv)
#Coefficients du modèle retenus par la pénalisation LASSO
coef(lasso.r.cv, s = "lambda.1se")
Résultats :
x.mat.test <- model.matrix(Energy ~ . - 1, data = test_reg)
lasso.r.cv.pred <- predict(lasso.r.cv, s = "lambda.1se", newx = x.mat.test) #valeur predite sur le jeu de train
Tab_score_reg <- Compute_Error(lasso.r.cv.pred,quali = FALSE, name_model = "Lasso penalisation");Tab_score_reg
Tab_F1_reg <- Compute_F1(lasso.r.cv.pred,quali=FALSE,name_model = "Reg LASSO")
Rappel : La pénalisation LASSO permet de diminuer le nombre de variable utilisé tout en dégradant très peu l'information. Dans ce modèle une contrainte est ajouté sur les paramètres. La pénalisation LASSO a donc une RMSE supérieur au modèle linéaire car le modèle linéaire minimise uniquement les moindres carrés. Nous avons donc un moins bon modèle mais moins de variables sont nécessaires pour le construire.
c. Régression linéaire avec sélection de variable par le critère Akaike Information Criterion (AIC) :¶
Mise en place du modèle :
#Test sur le modèle sans les intéractions
glm.r <- glm(Energy~ ., data=train_reg)
glm.r.step <- step(glm.r, direction = "backward", trace=0)
anova(glm.r.step, test='F')
Résultats :
glm.r.step.pred <- predict(glm.r.step, newdata=test_reg)
Tab_score_reg <- Compute_Error(glm.r.step.pred,quali = FALSE, name_model = "Regression Linéaire Selection Variable");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.r.step.pred,quali=FALSE,name_model = "Reg AIC")
Mise en place du modèle quadratique :
glm.2.r <- glm(Energy~ .^2, data=train_reg)
# du critère d'Akaïke par méthode descendante
glm.2.r.step <- step(glm.2.r, direction = "backward", trace=0)
#On réalise un test de fischer sur le résultat pour confirmer que toutes les variables retenues sont bien significatives
anova(glm.2.r.step, test='F')
Résultats :
glm.2.r.step.pred <- predict(glm.2.r.step, newdata=test_reg)
Tab_score_reg <- Compute_Error(glm.2.r.step.pred,quali = FALSE, name_model = "Regression Linéaire Quadratique Selection Variable");Tab_score_reg
Tab_F1_reg <- Compute_F1(glm.2.r.step.pred,quali=FALSE,name_model = "Reg AIC quadra")
Commentaires : Deux regressions linéaires avec selection de variables ont été réalisées : le modèle initial et un modèle quadratique. La selection de variables (ici méthode descendante) permet de garder les variables les plus intéressantes au sens du critère AIC.
En comparant les différents scores, la regression quadratique avec la selection de variables est très légèrement meilleure.
d. Arbre de décision optimal¶
Mise en place du modèle :
options(repr.plot.width = 15, repr.plot.height = 10)
tree.r=rpart(Energy~.,data=train_reg,control=rpart.control(cp=0.001))
xmat=xpred.rpart(tree.r, xval = 10)
#levels(datappr[,"Survived"]) <- c(1,2)
xerr <- (xmat-as.numeric(y_train_reg))^2
CVerr=apply(xerr,2,sum)
#CVerr # CP erreur
cp_opti = as.numeric(attributes(which.min(CVerr))$names);
cat('cp opti :', cp_opti)
tree.opti.r=rpart(Energy~.,data=train_reg,control=rpart.control(cp=cp_opti))
options(repr.plot.width = 15, repr.plot.height = 10)
plot(tree.opti.r)
text(tree.opti.r)
Commentaire : On remarque que, comme pour le modèle de classification, la variable Overall.height a un grand impact sur le modèle.
Résultats :
tree.opti.r.pred <- predict(tree.opti.r, newdata = test_reg)
Tab_score_reg <- Compute_Error(tree.opti.r.pred,quali = FALSE, name_model = "Tree Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(tree.opti.r.pred,quali=FALSE,name_model = "TreeOpti Reg")
Commentaire :
Tout comme précedemment en régression, nous avons une amélioration significative des résultats avec un arbre de décision (méthode non linéaire).
e. Random Forest¶
Mise en place du modèle :
set.seed(69)
rf.r <-randomForest::randomForest(x = x_train_reg, y = y_train_reg,
xtest = x_test_reg, ytest = y_test_reg,
importance = TRUE)
Résultats :
rf.r.pred <- rf.r$test$predicted
Tab_score_reg <- Compute_Error(rf.r.pred,quali = FALSE, name_model = "Random Forest");Tab_score_reg
Tab_F1_reg <- Compute_F1(rf.r.pred,quali=FALSE,name_model = "Reg RF")
Commentaire :
Les résultats obtenus avec Random Forest sont similaires à l'arbre de décision optimisé. Cependant, la RMSE est meilleur pour Random Forest.
Tout comme précemment, nous allons tenter d'améliorer le modèle Random Forest en optimisant le paramètre mtry correspondant au nombre de variables testées à chaque division.
Mise en place du modèle optimisé :
set.seed(69)
cvControl <- trainControl(method = "cv", number = 10)
mtryTrials <- train(train_reg, y_train_reg, method = "rf", tuneLength = 7,
trControl = cvControl, trace = FALSE)
set.seed(69)
defaultW <- getOption("warn")
options(warn = -1)
rf.opti.r<-randomForest::randomForest(x = x_train_reg, y = y_train_reg,
xtest = x_test_reg, ytest = y_test_reg,
mtry=mtryTrials$bestTune$mtry,importance = TRUE)
options(warn = defaultW)
Résultats :
rf.opti.r.pred <- rf.opti.r$test$predicted
Tab_score_reg <- Compute_Error(rf.opti.r.pred,quali = FALSE, name_model = "Random Forest Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(rf.opti.r.pred,quali=FALSE,name_model = "Reg RF Opti")
Commentaire :
Le modèle Random Forest retourne jusqu'ici les meilleurs résultats en terme de RMSE.
Cependant, par rapport aux autres métriques, les modèles Random Forest et d'arbre de décision sont équivalents.
f. Boosting¶
Mise en place du modèle :
gbm.r=gbm(Energy~., data=train_reg,distribution="gaussian",n.trees=2000, cv.folds=10,
n.minobsinnode = 5,shrinkage=0.05,verbose=FALSE)
plot(gbm.r$cv.error)
# nombre optimal d'itérations par valiation croisée
best.iter=gbm.perf(gbm.r,method="cv")
gbm.r.pred=predict(gbm.r,newdata=test_reg,n.trees=best.iter)
Tab_score_reg <- Compute_Error(gbm.r.pred,quali = FALSE, name_model = "GBM Naive");Tab_score_reg
Tab_F1_reg <- Compute_F1(gbm.r.pred,quali=FALSE,name_model = "Reg GBM")
Mise en place du modèle optimisé:
cvControlRandom = trainControl(method = "cv", number = 10, search='random')
gbm.tune.r = train(x_train_reg, y_train_reg, method='gbm', tuneLength=4, trControl=cvControlRandom, verbose=FALSE)
Résultats :
gbm.tune.r.pred = predict(gbm.tune.r, newdata = x_test_reg)
Tab_score_reg <- Compute_Error(gbm.tune.r.pred, quali = FALSE,
name_model = "GBM Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(gbm.tune.r.pred,quali=FALSE,name_model = "Reg GBM Opti")
Commentaire :
Tout comme avec la classification, la méthode de Gradient Boosting ne semble pas apporter de réelles améliorations par rapport à Random Forest, même avec des paramètres optimisés.
g. SVM¶
Mise en place du modèle optimisé:
svm.tune.r = tune.svm(Energy~.,data=train_reg,
cost=c(1,1.25,1.5,1.75,2),gamma = seq(0.1, 2, by = 0.2))
Commentaires :
Etant limité par les performances de nos machines, nous ne pouvons pas tester plusieurs noyaux pour optimiser le modèle SVM en régression.
Résultats :
svm.opti.r.pred=predict(svm.tune.r$best.model,newdata=test_reg)
Tab_score_reg <- Compute_Error(svm.opti.r.pred,quali = FALSE, name_model = "SVM Opti");Tab_score_reg
Tab_F1_reg <- Compute_F1(svm.opti.r.pred,quali=FALSE,name_model = "Reg SVM Opti")
Commentaire :
Nous n'obtenons pas de meilleur résultat avec SVM, de plus cette méthode est très couteuse.
Bilan - Régression¶
cat("Bilan des scores en passant par la Regression : \n " )
Tab_score_reg
Pour l'interprétation, se référer à la partie Conclusion.
D. Conclusion¶
cat("Bilan des scores - Classification : \n " )
Tab_score_class
cat("Bilan des scores - Regression : \n " )
Tab_score_reg
Pour rappel, nous pouvons considérer un modèle comme étant performant si il a une forte accuracy, un faible MAE et un MacroF1 élevé. Dans le cas d'une régression, on considère la RMSE comme étant un bon indicateur.
En classification, le meilleur modèle obtenu est celui de Random Forest : il est gagnant sur toutes les métriques.
En régression, nous observons une équivalence des modèles Random Forest et Arbres de décision : Random Forest a une RMSE plus faible mais une Macro F1 moins bonne que celle de l'arbre de décision.
A première vue l'accuracy ne semble pas élevé, mais nous pouvons nuancer cette constatation par le fait que la classification se fait selon 7 classes. Aléatoirement, nous avons 1 chance sur 7 de choisir la bonne classe soit une "accuracy" aléatoire de 0.15. Ici notre modèle est donc nettement supérieure à cette probabilité.
De plus, les erreurs de classification sont souvent réalisées dans les classes voisines (cf MAE).
En prenant en compte tous ces éléments nos modèles semblent relativement pertinents.
Analyse des scores F1¶
Les graphiques suivants exposent les scores F1 de chaque classe énergétique pour chaque modèle implémenté.
l = nrow(Tab_F1_class)
c = ncol(Tab_F1_class)
gras = 3
plot(Tab_F1_class$Modele,rep(-1,l),type='l', ylim=c(0,1),
xlab="Modele", ylab= "F1 Score", main = "Score F1 pour les différentes classes énergétiques pour les modèles de classification")
grid()
points(Tab_F1_class$Modele,Tab_F1_class$A, col=1, lwd=gras)
points(Tab_F1_class$Modele,Tab_F1_class$B, col=2, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$C, col=3, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$D, col=4, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$E, col=5, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$F, col=6, lwd = gras)
points(Tab_F1_class$Modele,Tab_F1_class$G, col=7, lwd = gras)
legend(x='bottomright', legend=c("A","B","C","D","E","F","G"), pch='o', col=1:(c-1))
l = nrow(Tab_F1_reg)
c = ncol(Tab_F1_reg)
gras = 3
plot(Tab_F1_reg$Modele,rep(-1,l),type='l', ylim=c(0,1),
xlab="Modele", ylab= "F1 Score", main = "Score F1 pour les différentes classes énergétiques pour les modèles de régression")
grid()
points(Tab_F1_reg$Modele,Tab_F1_reg$A, col=1, lwd=gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$B, col=2, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$C, col=3, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$D, col=4, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$E, col=5, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$F, col=6, lwd = gras)
points(Tab_F1_reg$Modele,Tab_F1_reg$G, col=7, lwd = gras)
legend(x='bottomright', legend=c("A","B","C","D","E","F","G"), pch='o', col=1:(c-1))
En classification, ce sont les classes A et G qui sont les mieux prédites pour tous les modèles. Ce résultat est cohérent car ce sont des classes extrêmes avec de larges intervalles (A $\in [0,30[$ et G>75). En régression, cette constatation est également vraie hormis pour les modèles linéaires.
En classification, la classe B a du mal à être prédite. Ce résultat est aussi chohérent car l'intervalle de la classe B est inférieur aux autres (B $\in [30,35[$). Le risque de se tromper est donc augmenté. En régression, on ne distingue pas de classe majoritairement discriminée.
En mettant en parallèle les deux graphiques : il est interessant de constater pour Random Forest que les valeurs de macroF1 sont plus élevées en regression que en classification. Par ailleurs le delta en regression est moins important que celui en classification. On peut donc considérer le modèle de regression comme meilleur.
Une piste d'ouverture de notre étude consisterait à implémenter un algorithme permettant d'abord de prédire le type de consommation (Faible ou Forte) correspondant aux classes (A,B,C ou D,E,F,G) et ensuite d'appliquer un algo de classification sur ces données.
Lien entre les modèles et l'analyse statistique initiale :¶
Les résultats obtenues avec les différents modèles sont cohérents avec notre analyse exploratoire des données.
- Répartition des classes
On voit tout d'abord qus les classes B,C,D sont très peu représentées dans le jeu de données initial et que la classe A domine. Cela a un impact direct sur la classification, les individus auront tendance à être prédits dans la classe A, même si ils sont dans la classe B ou la classe C. Nos modèles apprennent sur moins de données provenant de ses classes, tout cela combiné à leur situation de classes intermédiaires (donc plus difficiles à prédire car non extrêmes) entraine une réelle difficulté à vraiment segmenter les classes B et C, comme le montrent nos graphiques représentant les métriques F1 selon les classes.
Une solution pourrait être de choisir des jeux d'apprentissages avec plus d'individus de ces classes mals représentées.
- Variables retenues dans les modèles
Ensuite on voit que dans nos modèles pénalisés les variables retenues sont, dans le cas de LASSO avec lambda=1SE, Overall.height, Glazing Area et Wall Area. Cela s'explique car les variables Wall.Area et Glazing.Area sont les variables les moins corrélées avec les autres dans notre matrice des corrélations. On ne peut donc pas s'en passer car leur information n'est liée à aucune autre variable et ne peut donc se retrouver ailleurs.
On observe aussi que la première variable servant à départager le premier noeud des arbres est overall.height. Dès les premiers graphiques montrant la hauteur selon la classe énergétique, on observait déjà que la variable Overall.height jouait un rôle important. En effet, les classes A,B et C sont toutes de hauteur 3,5m contre 7m pour les classes D,E,F et G. La deuxième variable intervenant dans les arbres pour la séparation est Glazing Area, variable portant à elle seule l'axe 2 de l'ACP (et expliquant 25% de la variance). Cette variable a donc évidemment un rôle important dans les modèles.
Enfin, le clustering ne mettait pas en évidence 7 clusters, ce qui laissait présager que la classification ne serait pas évidente. On remarque par ailleurs que la méthode du coude nous suggère 2 clusters (correspondant aux faibles et fortes consommations) et c'est réellement sur ce point que nous sommes bons : prédire une forte ou faible consommation (cf partie Bonus).
Enfin, les méthodes d'apprentissage non linéaires sont plus efficaces ici, même en regression.
Bonus - classification hierarchique¶
Nous avons vu de part le clustering et l'ACP, deux classes se distinguaient bien (ABC vs DEFG) nous allons donc mettre en place une classification hiérarchique :
1. Mise en place d'un modèle de prédiction binaire : Bonne ou mauvaise consommation ?¶
#On fait une nouvelle data_frame relabellisée :
x_bin = x_fin
Conso = ifelse(x_fin$Energy<=45,"Faible","Forte")
x_bin = cbind(x_bin,Conso)
x_bin$Conso <- factor(x_bin$Conso)
head(x_bin)
train_bin <- x_bin[train_ind, ]
test_bin <- x_bin[-train_ind, ]
train_bin <- dplyr::select(train_bin, Relative.compactness,Surface.area, Wall.area,Overall.height,orientation, Glazing.area,Glazing.area.distr,Conso)
test_bin <- dplyr::select(test_bin, Relative.compactness, Surface.area,Wall.area,Overall.height,orientation, Glazing.area,Glazing.area.distr,Conso)
prems <- glm(Conso ~., data=train_bin, family = binomial(link = "logit"))
res.fit <- predict(prems, newdata = test_bin, type="response")
t1=table(res.fit>0.5,test_bin$Conso);t1
score = sum(diag(t1)/sum(t1))
print(score)
2. Mise en place des modèles de prédiction pour les consos fortes et faibles¶
Séparation des données¶
train_bin <- x_bin[train_ind, ]
train_bin <- dplyr::select(train_bin,-Energy.efficiency)
test_bin <- x_bin[-train_ind, ]
test_bin <- dplyr::select(test_bin,-Energy.efficiency)
#juste : on sépare les forts et les faibles reels pour entrainer nos modèles
train_bin_Forte <-train_bin[which(train_bin$Conso=='Forte'),]
train_bin_Faible <-train_bin[which(train_bin$Conso=='Faible'),]
#vraies valeurs correspondant à nos predictions de bonne ou mauvaise classe
test_bin_Forte <- test_bin[which(res.fit > 0.5),]
test_bin_Faible <- test_bin[which(res.fit < 0.5),]
#On enlève la variable Conso :
train_bin_Forte <- dplyr::select(train_bin_Forte,-Conso)
train_bin_Faible <- dplyr::select(train_bin_Faible,-Conso)
test_bin_Forte <- dplyr::select(test_bin_Forte, -Conso)
test_bin_Faible <- dplyr::select(test_bin_Faible,-Conso)
x_train_bin_Forte <- dplyr::select(train_bin_Forte, -Energy)
y_train_bin_Forte <- train_bin_Forte$Energy
x_test_bin_Forte <- dplyr::select(test_bin_Forte, -Energy)
y_test_bin_Forte <- test_bin_Forte$Energy
x_train_bin_Faible <- dplyr::select(train_bin_Faible, -Energy)
y_train_bin_Faible <- train_bin_Faible$Energy
x_test_bin_Faible <- dplyr::select(test_bin_Faible, -Energy)
y_test_bin_Faible <- test_bin_Faible$Energy
Modele de prediction : fortes consommations¶
rf.h.forte <- randomForest(x = x_train_bin_Forte, y = y_train_bin_Forte ,
xtest = x_test_bin_Forte, ytest = y_test_bin_Forte,
ntree = 500,do.trace = 50, importance = TRUE)
rf.h.forte.pred <- rf.h.forte$test$predicted
class.forte.pred<- classify(rf.h.forte.pred)
y.forte.class <- classify(y_test_bin_Forte)
table.forte = table(pred.reg = class.forte.pred , observations = y.forte.class); table.forte
Modele de prediction : faibles consommations¶
rf.h.faible <- randomForest(x = x_train_bin_Faible, y = y_train_bin_Faible ,
xtest = x_test_bin_Faible, ytest = y_test_bin_Faible,
ntree = 500,do.trace = 50, importance = TRUE)
rf.h.faible.pred <- rf.h.faible$test$predicted
class.faible.pred<- classify(rf.h.faible.pred)
y.faible.class <- classify(y_test_bin_Faible)
table.faible = table(pred.reg = class.faible.pred , observations = y.faible.class); table.faible
3. Mise en commun des modèles¶
pred.hierarc = c(class.faible.pred,class.forte.pred)
real.hierarc = c(y.faible.class,y.forte.class )
table.hierarc = table(pred.reg = pred.hierarc, observations = real.hierarc); table.hierarc
4. Resultats¶
cat("Score : ",sum(diag(table.hierarc))/sum(table.hierarc), "\n" )
cat("Mean Absolute Error : ", MAE(table.hierarc), "\n")
cat("Macro F1 : ", mean(macro_F1(table.hierarc)) )
Bilan : Les scores obtenus sont bons, juste en dessous de notre meilleur prédicteur des parties précédentes qui était le random Forest.
Premièrement, le modèle de régression logistique binaire pour prédire entre les consommations d'energie fortes et faibles est très bon (précision de 0.96).
En regardant la matrice de confusion, on voit que cette méthode de classification deux modèles qui apprennent spécifiquement sur des données de consommation faible ou fortes, n'apporte pas un tel gain car on se trompe encore beaucoup pour les classes intermédiaires (B,C, D et E). Plus particulièrment, on voit que la prédiction pour un individu se situant réellement dans la classe B est toujours compliquée, on va se tromper 2 fois sur 3 de prédiction en prédisant soit la classe A ou C.
On peut faire le lien avec ce que nous avons dit précedemment dans la conclusion, ces erreurs de prédictions sont dues aux intervalles choisies pour appartenir aux classes.
On peut finalement se demander si un tel modèle ne serait pas plus performant si on avait plus de données. Il se pourrait en effet que les modèles basés sur les consommations fortes et faibles soient un peu plus performants avec plus d'apprentissage.